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ABSTRACT 


The  problem  of  generating  random  Toeplitz  correlation  matrices  is  considered. 
Several  methods  are  proposed,  of  which  the  most  promising,  as  determined  by  both 
computational  complexity  and  spectral  randomness,  seems  to  be  that  based  on  the 
characteristic  function  of  random  discrete  probability  measure.  A  number  of  inter¬ 
esting  theoretical  issues  are  recorded  for  further  investigation.  This  report  follows 
an  earlier  one  devoted  to  general  random  correlation  matrices.  In  both  ca^es  the 
intent  is  to  use  such  matrices  to  simulate  random  data  for  the  testing  of  certain 
group-theoretic  signal  processing  algorithms. 
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1.  INTRODUCTION 


The  following  is  intended  as  a  sequel  to  our  earlier  report  on  random  correlation  matrices  [11]. 
All  this  work  is  ultimately  motivated  by  eventual  application  to  the  testing  and  comparison  of 
certain  signal  processing  algorithms,  specifically  those  defined  by  group  transforms  and  filters,  and 
intended  for  such  second  order  tasks  as  data  compression  and  decorrelation,  and  Wiener  filtering. 
The  earlier  report,  particularly  its  summary,  contains  further  details  on  these  applications. 

The  basic  idea  is  that  a  random  correlation  matrix  (RCM)  is  to  be  interpreted  as  the  covariance 
matrix  of  a  random,  or  ‘average,’  signal.  The  latter  is  conceived  as  a  random  vector,  whose 
components  all  have  an  equal  variance  (here  taken  as  unity).  This  vector  arises,  in  turn,  from 
sampling  and  quantizing  a  continuous  signal.  If  the  latter  can  be  reasonably  modeled  as  a  stationary 
stochastic  process,  then  the  resulting  covariance  matrix  of  the  samples  will  have  a  Toeplitz  structure, 
that  is,  it  will  be  Hermitian  and  have  constant  diagonals.  Hence  it  can  be  fully  specified  by,  for 
example,  its  first  row. 

Our  view  is  that  a  Toeplitz  random  correlation  matrix  (TRCM)  represents  an  intermediate 
degree  of  statistical  structure  between  general  random  signals  as  described  by  ROMs,  and  signals 
describable  as  outputs  of  very  special  models,  such  as  linear  stochastic  difference  equations.  The 
latter  result  in  autoregressive  signal  models  of  some  finite  order  m,  denoted  AR{Tn).  In  our  numer¬ 
ical  studies  of  group  -  theoretic  signal  processing  algorithms,  we  restrict  to  the  cases  m  <  2.  There 
is  an  interesting  cluster  of  theoretical  questions  here,  pertaining  to  the  density  of  such  signals  and 
their  matrices,  which  will  be  listed  later  in  Section  4,  along  with  other  unresolved  issues. 

The  specific  problem  which  we  address  below  is  that  of  generating  TRCMs.  As  in  the  earlier 
report,  we  distinguish  two  aspects  of  this  problem:  analytical  and  computational.  The  latter  refers 
to  computer  software  whose  implementation  results  in  an  output  stream  of  pseudorandom  TRCMs. 
We  will  also  have  a  few  remarks  to  make  about  some  spectral  distributions  associated  with  our 
TRCMs,  but  not  to  the  extent  of  those  in  [llj;  most  of  the  natural  questions  in  this  context  are 
left  to  Section  4. 

The  general  topic  of  Toeplitz  matrices  and  operators  is  extensively  developed  by  now,  and  in 
no  way  encompassed  by  a  few  sources.  Some  general  references  are  the  books  [10,  12]  and  articles 
[9,  19],  but  these  already  provide  more  background  than  is  really  necessary  for  present  purposes. 

The  subject  of  random  correlation  matrices  as  presented  in  this  and  the  preceding  report  [11], 
and  in  several  references  listed  there,  is  but  one  aspect  of  the  much  larger  area  of  random  matrices. 
As  we  see  it,  this  area  encompasses,  in  addition  to  the  present  topic,  the  asymptotics  of  the  spectral 
distributions  of  large  random  matrices,  given  distributional  information  about  the  matrix  entries, 
and  the  problem  of  sampling  from  various  compact  matrix  groups  (symmetric  groups,  orthogonal 
groups,  etc.).  Collectively,  random  matrix  theory  has  a  vast  array  of  applications  in  the  sciences 
(nuclear  physics,  biology),  communications  engineering  (coding,  encryption),  statistics  (simulation, 
projection  pursuit,  algorithm  tests),  and  mathematics  (random  walks,  spectral  asymptotics,  ergodic 
theory).  The  recent  conference  proceedings  [16]  gives  an  overview  of  many  of  these  applications, 
along  with  considerable  theory.  The  books  [3,  18]  give,  respectively,  a  physical  and  a  mathematics 
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perspective  on  random  matrices,  with  emphcisis  on  eigenvalue  distributions.  For  the  early  work, 
emanating  from  the  need  to  statistically  model  energy  levels  of  complex  nuclei,  dating  at  lea^st  from 
E.  Wigner  (1957),  see  his  survey  article  [21].  Finally,  there  are  two  Russian  books  by  V.  Girko 
(Kiev,  1975  and  1980)  on  random  matrices  and  random  determinants,  respectively,  which  have  not 
been  available  to  the  author. 
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2.  THE  METHODS 


In  this  section  we  will  discuss  five  possible  methods  for  defining  TRCMs.  The  first,  an  attempt 
at  extending  the  random  Gram  method  of  [11]  for  generating  ROMs,  is  rejected  as  being  infeasible. 
The  second  method,  derived  from  the  well-known  structure  of  moving  average  processes,  is  simple 
•and  efficient,  and  was  our  initial  choice.  The  third,  which  builds  up  a  matrix  from  the  character¬ 
istic  function  of  a  random  probability  distribution,  is  the  most  interesting  theoretically,  it  can  be 
reasonably  efficient,  and  can  give  good  spectral  distribution  properties,  at  least  for  small  matrix 
size.  The  remaining  two  have  various  deficiencies  of  both  a  theoretical  and  practical  nature  that 
probably  render  them  unsuitable  for  automatic  use.  Yet  each  of  them  involves  interesting  open 
questions,  and  thus  it  is  felt  that  they  are  worthy  of  brief  mention  in  the  present  context. 


Before  detailing  the  five  methods,  let  us  be  specific  about  what  we  are  trying  to  do.  For  a 
fixed  positive  integer  N  we  want  a  procedure  for  producing  a  stream  of  x  matrices  of  the  form 


1  a  b  . . .  c 

a  1  a 

b  a  1 

:  ...  a 

c  .  a  1 


(2.1) 


each  of  which  is,  in  addition,  positive  definite.  The  N-  1  free  entries  b.,  . . .  are  to  be  deterministic 
functions  of  some  random  variables.  In  the  real  symmetric  case,  of  course,  we  will  have  a  =  a, 
6  =  6,  . . . ,  etc.  For  computer  implementation  we  will  expect  to  be  able  to  reduce  all  randomness 
to  pseudorandom  samples  from  the  univariate  uniform  or  normal  distributions.  Any  matrix  of  the 
form  (2.1)  constructed  by  some  such  procedure  is  by  definition  a  TRCM. 


Method  1  This  is  based  on  the  procedure  extensively  discussed  in  [11).  That  procedure 
yields  what  were  called  there  random  Gram  matrices.  Namely,  for  fixed  A^,  we  select  N  random 
vectors  ri, . . . ,  ryv  independent  and  uniformly  distributed  on  the  unit  sphere  in  the  real  space  of  N 
(or  more)  dimensions,  and  form  the  Gram  matrix 

=  [(t;j,Vj)].  (2.2) 

General  theory  forces  A  to  be  almost  surely  positive  definite  and  hence  A  is  a  correlation  matrix, 
since 


Extension  of  this  method  to  the  present  Toeplitz  case  then  reduces  to  the  question  of  whether 
this  further  symmetry  can  be  attained  by  some  geometrical  constraints  on  the  mutual  position  of 
the  vectors  Uj.  For  example,  in  the  N  =  3  case,  we  would  choose  V\,  V3  independent,  and  then  V2  at 
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random  from  the  unit  sphere  in  the  plane  (t;i  —1^3)"^.  This  would  guarantee  that 
and  hence  that  the  Gram  matrix  A  in  Equation  (2.2)  is  Toeplitz.  With  a  little  more  effort  an 
analogous  procedure  could  be  developed  for  the  =  4  case,  but  we  have  not  been  able  to  obtain 
such  an  algorithm  for  general  N. 


Method  2  This  is  based  on  the  well-known  formula  for  the  autocorrelation  function  of 
a  (stationary)  moving  average  process.  If  {x^}  is  such  a  process,  of  order  p,  then  there  are  constants 
cq  =  1,  Cl, ...,  Cp,  such  that  xt  =  wt  CiWt-i  -h  •  •  •  +  CpWt-p^  where  is  a  white  noise  process. 
It  follows  that 


cov{xt+k,Xt)  =  (2.3) 

j=0 

for  |fc|  <  p,  and  0  otherwise.  Here  is  the  variance  of  the  noise  terms  Wt]  we  may  as  well  take 
(7  =  1  for  present  purposes.  The  function  of  k  defined  by  Equation  (2.3)  is  the  autocorrelation 

function  of  the  process  {x^},  up  to  normalization  by  the  factor  (1  -h  Cj  H - -(-  Cp),  and,  as  such,  is 

in  particular  a  positive  definite  function  on  the  group  Z_  of  integers. 

Now,  by  definition,  this  property  of  a  function  (j)  defined  on  Z_  is  exactly  that  required  for  the 
matrix  A  defined  by 


A  =  \4>{i  -  j)],  (2.4) 

to  be  a  Toeplitz  covariance  matrix.  Hereafter  we  will  enforce  the  normalization  <^(0)  =  1,  so  that 
A  is  actually  a  correlation  matrix.  So,  at  this  point  we  may  say  that  a  normalized  positive  definite 
function  on  Z_  will  give  rise  to  a  Toeplitz  correlation  matrix  of  any  order  through  the  prescription 
of  Equation  (2.4),  and  that,  for  a  fixed  order,  a  simple  way  to  obtain  such  functions  is  through  the 
covariance  formula  (2.3). 

With  this  background  we  can  now  state  a  simple  algorithm  for  generating  real  TRCMs: 

•  Select  matrix  dimension  N\ 

•  Generate  a  vector  (xi, . . . ,  xyv)  from  the  uniform  distribution  on  the  unit  sphere 
in  \ 

•  Define  a  function  <{>  on  by  0(/c)  =  XjXj^k  with  the  convention  that 

Xj  =  0,  j  >  N; 

•  Define  the  associated  correlation  matrix  C  by  C  =  [<t>{\i  —  j|)]. 

We  note  that  procedures  for  carrying  out  Step  2  above  are  discussed  in  [11,  Section  5.1)  and, 
in  greater  depth,  in  [6,  Section  5.4).  Some  spectral  properties  of  this  method  are  given  in  the  next 
section. 
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Method  3  We  continue  with  the  theme  of  using  positive  definite  functions  to  generate 
Toeplitz  correlation  matrices  according  to  Equation  (2.4).  Hence  we  will  permit  an  interest  in 
generating  complex  Hermitian  Toeplitz  correlation  matrices,  as  displayed  in  Equation  (2,1).  (If 
real  Toeplitz  correlation  matrices  are  required,  then  the  real  parts  of  the  eventual  complex  matrices 
may  be  taken;  doing  so  will  not  alter  the  spectrum.)  For  this  it  is  sufficient  to  have  a  procedure  for 
generating  positive  definite  functions  on  the  group  R  of  real  numbers.  Of  course,  there  is  a  slight 
catse  of  ‘overkill’  in  this  approach  since  such  a  function  (f)  will  produce  infinitely  many  matrices  of 
the  desired  sort,  and  of  any  dimension  N,  by  the  formula 


C  =  [</.(<,■  -  <,)], 


(2.5) 


where  {^i}  C  i?,  i,  j  =  1, . . . ,  A^,  and  we  would  want  the  work  required  to  specify  <p  to  be  somehow 
scaled  to  the  size  of  N. 

We  make  use  of  Bochner’s  theorem  which  states  that  each  continuous  normalized  positive 
definite  function  is  the  characteristic  function  of  a  probability  distribution  on  R.  That  is 

(^(x)=  [  €^^^dP{t),  (2.6) 

Jr 

where  P{^)  is  a  distribution  function  on  R^  or,  equivalently,  a  probability  measure  on  (the  Borel 
subsets  of)  R.  Hence,  via  the  correspondence  P  (l>  we  see  that  a  method  for  generating 

TRCMs  will  follow  from  any  method  that  randomly  selects  a  probability  measure  P  on  R. 

Such  a  method  would  seem  to  involve  us  with  the  concept  of  a  probability  measure  over  the 
(enormous!)  space  Il{R)  of  all  probability  measures  on  R.  In  fact,  such  measures  have  been 
introduced  in  the  context  on  nonparametric  Bayesian  statistical  analysis  by  T.  Ferguson  in  1973- 
74,  and  termed  Dirichlet  processes.  (Actually,  the  concept  of  a  random  probability  measure  with 
fixed,  bounded  support,  e.g.,  [0,  l],  identified  with  the  corresponding  distribution  function,  goes 
back  at  least  10  years  earlier  [7].)  To  engage  in  a  substantial  discussion  of  these  would  constitute 
a  further  case  of  ‘overkill’  in  the  present  context.  So  let  us  just  say  that  Dirichlet  processes  can 
be  parameterized  by  the  set  of  finite  measures  on  R^  that  is,  measures  of  the  form  /3  =  ip^t  > 
0,p£ll{R),  If  Dp  is  the  Dirichlet  process  corresponding  to  such  a  /?,  then  for  any  partition  R  = 
BiU. .  .UjBfc,  the  random  vector  [Dp{Bi)^ . . . ,  Dp{Bk)]  has  the  Dirichlet  distribution  with  parameters 
(/3(jBi),.  . .  ,/3(jBfc)).  In  particular,  we  can  say  that  E{D0)  =  p{=  t~^P)  .  A  basic  survey  article  is  [8] 
with  other  early  work  by  C.  Antoniak,  D.  Blackwell,  and  J.  MacQueen.  There  have  been  sporadic 
developments  since  then  associated  with  such  authors  as  S.  Dalai,  J.  Sethuramen,  R.  Tiwari,  inter 
alia.  Some  of  this  work  involves  the  concept  of  an  exchangable  sequence  of  random  variables. 

One  key  result  due  independently  to  several  of  the  early  authors  is  that  a  p£ll{R)  chosen 
at  random  from  a  Dirichlet  process  is  almost  surely  discrete.  A  particularly  nice  constructive 
specification  of  a  Dirichlet  process  was  given  in  [2].  It  is  based  on  abstract  ‘Polya  urn  scheme’ 
(drawing  balls  of  various  colors  from  an  urn,  with  a  prior  distribution  on  the  colors;  infinitely  many 
colors  allowed),  and  yields  approximating  measures  of  the  form 
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(2.7) 


with  =  0{R)  +  N ^  and  (5(i)  =  the  point  measure  concentrated  at  tsR-  It  is  clear  that  such  mea¬ 
sures  could  be  rather  easily  simulated  on  a  computer.  From  there,  the  resulting  is  substituted 
for  P  in  Equation  (2.6),  thereby  defining  the  function  0  which,  in  turn,  is  used  in  Equation  (2.5) 
to  define  a  Toeplitz  correlation  matrix  C  (with  some  particular  choice  of  points  U). 

We  have  not  carried  out  this  program  although  it  appears  to  be  a  very  interesting  project, 
particularly  because  of  the  likely  theoretical  realtionships  between  various  spectral  properties  of  the 
resulting  TRCMs  and  the  probabilistic  structure  of  the  random  measures  that  might  thereby 
be  suggested. 

Instead,  we  offer  the  following  ad  hoc  algorithm  for  generating  a  discrete  random  probability 
measure  P,  which  is  then  used  in  the  usual  way  to  define  a  Toeplitz  correlation  matrix: 

•  Select  matrix  dimension  N; 

•  Select  an  integer  L  >  N ^  and  a  positive  number  a; 

•  Generate  a  random  point  (pi, . . .  ,pl)  from  the  unit  simplex,  so  that 
Pi  >  0,  Epi  =  1; 

•  Generate  a  random  point  (5i,...,s/,)  from  the  centered  multivariate  normal 
distribution  with  L  x  L  covariance  matrix  = 

•  Define  P  =  pi(5(si)  -h  . . .  +  Pl6(si)  ; 

•  Define  (t>  from  Equation  (2.5),  and  finally  C  from  Equation  (2.4)  with  U  =  i 
there. 

The  choice  of  the  integer  L  in  the  second  step  above  is  important,  as  too  large  or  too  small  a 
value  causes  a  distinctly  nonrandom  distribution  of  eigenvalues.  From  our  numerical  work  reported 
in  the  next  section,  it  appears  that  L  =  3N  is  a  good  choice.  Note  that  L  >  is  necessary  for 
positive  definiteness  of  C.  The  choice  of  a  is  less  crucial;  a£:[5A^,  lOA^]  seems  to  be  effective. 

A  procedure  for  carrying  out  the  third  step  above  was  discussed  in  [11,  Section  4.1];  it  reduces 
to  taking  p,-  to  be  the  zth  spacing  in  a  sample  of  size  L  —  1  from  the  standard  uniform  f/[0,  l] 
distribution. 

Method  4  This  is  based  in  the  geometric  notion  of  orthogonal  projection  in  matrix 
space.  We  consider,  in  the  real  or  complex  case  respectively,  the  real  ordered  vector  space  of  real 
symmetric  or  complex  Hermitian  N  x  N  matrices.  Let  us  call  this  space  ATyv  and  denote  its  linear 
subspace  of  all  Toeplitz  matrices  of  the  form  (2.1)  by  T/v-  Equipped  with  the  usual  Hilbert  - 
Schmidt  inner  product 


{A,B)  =  tr{ABn, 


(2.8) 


Xj\^  becomes  a  real  Hilbert  space,  and  we  can  contemplate  the  operation  of  orthogonal  projection 
of  Xj\[  onto  the  subspace  .  Let  us  denote  this  operation  by 

A  ^  At,  AeXr,  .  (2.9) 

It  turns  out  to  be  easy  to  compute  this  projection  because  of  the  presence  of  an  obvious  orthogonal 
basis  for  Namely,  if  we  specify  a  Toeplitz  symmetric  or  Hermitian  matrix  by  giving  its  first 
row,  then  the  basis  is  given  by  the  standard  unit  vector  basis  for  in  the  real  case,  and  these 
vectors  plus  y/  —  1  times  these  vectors  in  the  complex  case  (omitting  (1,0,...,  0)^  in  this  case). 
Thus 


dim(7iv)  = 


N,  real  case 
2N  —  1,  complex  case. 


If  we  examine  the  nature  of  the  projection  operation  by  computing  the  Fourier  coefficients  of 
a  matrix  in  X}\[^  At  is  obtained  by  averaging  each  diagonal.  Thus,  for  instance,  if  =  4,  and 


then 


A  = 


a  e  h  I 
b  f  k 
c  9 
d 


At  = 


(X  (3  j  6 

(X  P  j 

OC  P 


X 


with  oc=  (a  6  c  +  d)/4,  P  =  (e  +  /  +  ^)/3,  7  =  {h-\-  k)/2,  and  8  —  L 

It  follows  that  if  we  start  with  a  correlation  matrix  CeXj^^  its  Toeplitz  projection  Ct  will 
again  be  a  correlation  matrix,  provided  that  it  has  no  negative  eigenvalues.  Unfortunately,  this 
need  not  be  the  case.  Table  2-1  below  exhibits  for  TV  =  4,  5,  correlation  matrices  in  JVat,  their 
Toeplitz  projections,  and  a  negative  eigenvalue  for  each  of  the  latter.  Thus  we  cannot  conclude  that 
a  correlation  matrix  necessarily  projects  onto  a  Toeplitz  correlation  matrix  (surprisingly,  however, 
we  were  unable  to  find  a  counterexample  for  TV  =  3). 

In  spite  of  this  somewhat  disappointing  discovery,  an  interesting  empirical  phenomenon  occurs 
in  this  context,  one  that  seems  worthy  of  further  investigation  and  explanation.  Namely,  the 
empirical  frequency  of  the  outcome  (correlation  matrix  with  indefinite  Toeplitz  projection)  tends 
rather  rapidly  to  0  with  increasing  values  of  TV.  Evidence  for  this  statement  is  summarized  in 
Table  2-2  below.  The  underlying  experiment  was  conducted  by  generating,  for  fixed  TV,  several 
hundred  ROMs  of  Gram  type,  C  =  TT*,  with  T  an  TV  x  (TV -hi)  matrix  with  row  vectors  distributed 
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TABLE  2-1 

Correlation  Matrices  with  Indefinite 

Toeplitz  Projections 

N  =  4 

C 

min  a  (Cj) 

1.0.697  -.709  -.911 

1.0  -.747-702 

1.0  .611 

1.0 

1.0. 187 -.705  -.911 

-.117 

N  =  5 

1.0  .252  -.168  -.493  -.904 
1.0  -.614  .583 -.239 
1.0  -.613 -.203 
1.0  .353 
1.0 

1.0 -.155  .206  -.366  -.904 

-.144 

independently  and  uniformly  over  the  sphere  in  The  resulting  matrices  were  then 

projected  onto  TV  and  their  eigenvalues  computed  via  IMSL  routine  EIGRS,  an  eigen-program 
designed  for  real  symmetric  matrices. 

These  findings  seem  rather  curious  and  deserving  of  further  study.  We  note  that  the  set  r(A^) 
of  real  N  x  N  correlation  matrices  is  a  compact  convex  set  of  dimension  N{N  —  l)/2.  The  fact  that 
this  is  rapidly  increasing  with  N  suggests  that  perhaps  we  do  not  have  a  sufficiently  large  sample 
size  for  JV  >  16  in  Table  2-2  to  be  able  to  conclude  much  about  the  ‘true’  probability  that  an  x 
RCM  has  a  Toeplitz  correlation  matrix  as  its  projection.  But  the  fact  that  these  projections  can 
occasionally  fail  to  be  positive  definite  means  that  the  following  simple  algorithm  for  generating 
TRCMs  cannot  be  guaranteed  failproof: 

•  Generate  a  RCM  C; 

•  Return 

For  this  reason  it  cannot  be  recommended  for  numerical  work. 


Method  5  As  discussed  at  some  length  in  [11],  one  natural  criterion  for  a  class  of 
RCMs  to  be  ‘truly  random’  is  that  the  associated  class  of  spectra  be  uniformly  distributed  over  the 
simplex  S]\[  =  {Aj  >  0  :  EA^  =  A^}  of  appropriate  dimension.  We  will  discuss  some  tests  of  such 
randomness  for  the  TRCMs  produced  by  Methods  2  and  3  in  the  next  section.  But  first  we  want 
to  take  note  of  a  recently  proposed  algorithm  which  might  permit  a  direct  solution  of  the  problem 
of  generating  TRCMs  with  random  spectrum. 
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TABLE  2-2 

Empirical  Failure  Rate  of  RCM  to 

Project  onto  TRCM 

No.  Trials 

No.  Failures 

Failure  Rate 

N  =  3 

500 

0 

.000 

4 

500 

11 

.022 

5 

500 

8 

.016 

8 

500 

5 

.010 

16 

500 

0 

.000 

24 

300 

0 

.000 

32 

200 

0 

.000 

Since  it  is  easy  to  sample  randomly  (i.e.,  uniformly)  from  the  simplex  the  issue  quickly 
reduces  to  the  possibility  of  constructing  a  Toeplitz  correlation  matrix  with  a  given  spectrum.  In 
general  (for  iV  >  5)  it  is  unknown  if  this  problem  always  has  a  solution;  this  is  the  ‘inverse  Toeplitz 
eigenproblem.’  Nevertheless,  a  recent  article  by  D.  Laurie  [13]  offers  a  numerical  procedure  which, 
in  spite  of  the  lack  of  a  convergence  proof,  is  claimed  to  solve  this  problem  in  practice,  in  fact,  to 
exhibit  quadratic  convergence. 

Let  {Ai}  C  :R  be  a  tentative  spectrum  of  a  real  symmetric  Toeplitz  matrix,  and  set  A  = 
diag[Ai, . . . ,  Ajvj.  If  T  is  such  a  matrix  then  there  must  be  an  orthogonal  matrix  Q  such  that 
T  =  QKQ*  or,  equivalently,  Q*TQ  =  A,  a  diagonal  matrix.  Laurie’s  algorithm  involves  a  succession 
of  alternate  orthogonal  diagonalizations  of  a  symmetric  Toeplitz  matrix,  followed  by  a  new  choice 
of  Toeplitz  matrix  so  that  its  diagonal,  under  the  current  orthogonal  basis,  agrees  with  A.  This 
latter  step  amounts  to  solving  an  x  linear  system,  arising  from  use  of  the  basis  vectors  for 
the  subspace  Tjv  discussed  earlier  under  Method  4.  Some  simplifications  are  possible  due  to  the 
reciprocal  structure  of  the  eigenvectors  of  a  Toeplitz  matrix.  Various  Toeplitz  starting  values  are 
suggested. 

In  spite  of  the  alleged  rapid  convergence  properties  of  this  procedure,  it  still  seems  like  a  lot 
of  work  to  undertake  to  obtain  TRCMs,  with  a  fair  potential  for  round-off  accumulation.  For 
this  reason  we  have  not  attempted  to  program  and  operate  this  method.  We  have  included  it  in 
this  discussion  because  it  does  present,  granting  convergence  and  round-off  control,  a  method  for 
producing  TRCMs  with  random  spectrum.  We  remark  that  the  inverse  Toeplitz  eigenproblem  has 
been  raised  in  an  earlier  RCM  context  by  Marsaglia  and  Olkin  [15]. 
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3.  NUMERICAL  EXPERIENCE 


In  this  section  we  summarize  some  computer  experiments  with  Methods  2,  3,  and  4  of  the 
previous  section.  There  is  clearly  any  number  of  questions  that  one  can  raise  concerning  the  distri¬ 
bution  of  various  spectral  functions  of  TRCMs  for  a  specific  probabilistic  generating  mechanism. 
However,  the  study  of  such  questions  was  never  our  goal  for  undertaking  this  w^ork.  Rather,  as 
discussed  in  the  introductory  section,  w^e  are  primarily  interested  in  a  faist  and  reliable  method  for 
generating  TRCMs  with  a  nicely  dispersed  spectrum,  that  is,  TRCMs  w'hose  spectrum  is  more  or 
less  uniformly  distributed  over  the  simplex  Sjs/. 

At  the  end  of  the  previous  section  we  noted  that  Laurie’s  algorithm  when  combined  w4th  a 
generator  for  sampling  from  the  uniform  distribution  over  Sjv  provides  the  most  direct  conceptual 
solution  of  this  problem,  yet  appears  not  to  be  ‘f2Lst,’  relative  to  the  other  methods.  So  what  we 
really  want  to  see  now  is  whether  some  of  the  other,  faster,  methods  can  be  said  to  yield  TRCMs 
wdth  a  random  spectrum. 

How  might  we  try  to  decide,  absent  a  definitive  theorem,  whether  a  given  class  of  TRCMs  has 
a  random  spectrum?  This  is  a  matter  that  was  discussed  at  some  length  in  our  earlier  report  [11]; 
how^ever,  the  central  issues  will  be  briefly  recalled  next. 

First,  we  say  that  an  NxN  random  correlation  matrix  A  has  a  random  spectrum  if  its  spectrum 
(t{A)^  viewed  as  a  point  in  the  simplex  Syv,  is  uniformly  distributed  over  S^.  Now,  one  way  to 
obtain  this  distribution  is  by  means  of  the  spacings  determined  by  a  sample  from  the  standard 
uniform  distribution  on  the  interval  [0,1].  Namely,  if  <  ^(TV^i)  are  the  order  statistics 

of  a  random  sample  drawn  from  U[0, 1],  and  we  put  U(o)  =  0,U(7V)  =  L  then  the  spacings  of  this 
sample  are  defined  by 


Si  =  i  ^  i  ^  A. 


The  spacings  are  positive  numbers  (with  probability  one),  and  sum  to  unity.  As  noted  in  [11],  it 
can  be  proved  that  the  vector  N’(5i,...,5;v)  is  uniformly  distributed  over  Spj. 

Now'  suppose  given,  for  some  fixed  A,  a  sample  of  size  n  of  N  x  N  TRCMs  generated  by 
some  particular  procedure;  for  example,  by  Method  2  or  3.  For  each  matrix  C  in  the  sample  we 
compute  its  spectrum  <t(C).  There  are  now  (at  legist)  two  general  types  of  tests  that  can  be  applied 
to  test  the  null  hypothesis  that  cr(C)  is  uniformly  distributed  over  Sp;.  The  first  type  consists  of 
computing  a  functional  f{a{C))  of  the  spectrum  whose  distribution,  under  the  null  hypothesis,  is 
known.  The  second  type  involves  a  transformation  T  :  Sp;  [0, 1];  we  evaluate  T{a{C))  obtaining 
a  finite  subset  of  [0, 1]  and  test  that  this  sample  is  from  [/[0, 1].  (Transformations  of  Sp;  onto  other 
simplicial  regions,  and  tests  against  uniformity  there  are  also  possible;  one  such  transformation  is 
simply  projection  on  the  first  A  —  1  coordinates.) 
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Examples  of  the  first  type  of  test  are  given  by  choosing  the  spectral  functions  /  as 


/(Ai,...,AaO  =  max(A,),  or 

=  min(Ai),  or  (3.1) 

=  \/(Ai -I- . . .  +  A^). 

These  functions  are,  respectively,  ||C||,  ||C“^||,  and  ||C||/r  (the  Frobenius  norm  of  the  matrix  C). 
An  example  of  the  second  type  of  test  is  given  by  choosing  the  transformation  T  as  T(Ai, . . . ,  \s)  = 
where 

x\  =  Ai/7V, 

X2  =  xi  +  \2INy  (3.2) 


xat-i  =  a:Ar_2  +  A(;v-1)/^- 

We  see  that  0  <  xi  <  X2  <  . . .  <  x^v-i  <  1,  and  that  the  associated  spacings  are  given  by 

Sk  =  1  <k<N. 

Under  the  null  hypothesis  Hq  of  a  random  (uniform)  spectrum,  each  of  the  random  variables 
defined  by  Equation  (3.1)  has  a  distribution  about  which  much  is  known;  see  the  discussion  in  (11, 
Sec.  4].  In  particular,  the  distribution  of  ||C"^||  is  available  in  closed  form: 

Pr(||C“'||  >  x)  =  (1  --  x)^-',  0  <  X  <  1,  (3.3) 

while  the  other  two  distributions  have  closed-form  asymptotic  expressions. 

If,  on  the  other  hand,  we  elect  to  utilize  the  transformation  T  :  (A)  — ♦  (x)  defined  by  Equation 
(3.2),  then  we  will  arrive  at  n  samples  (x)  which,  under  the  null  hypothesis,  are  each  drawn  from 
t/(0, 1].  Each  of  these  may  then  be  subjected  to  a  goodness-of-fit  test  for  uniformity,  of  which  there 
are  many.  Some  fraction  of  these  will  presumably  pass  such  a  test,  and  we  will  then  have  to  decide 
whether  this  fraction  is  large  enough  so  as  not  to  reject  Hq. 

A  third  alternative  is  to  view  each  sample  (x)  as  a  single  point  in  the  cube  (0, 1]^"^.  Under  Hq 
this  set  of  n  points  is  a  random  sample  from  {7(0, 1]^“^  and  a  multivariate  goodness-of-fit  test  for 
uniformity  could  be  applied.  This  approach  is  really  more  definitive  than  the  one  just  discussed  but, 
especially  for  large  A^,  would  require  much  larger  sample  sizes  than  what  are  presently  available. 

How  are  we  to  choose  among  this  multitude  of  plausible  tests?  Classical  statistical  theory 
suggests  that  we  consider  power.  But  this  too  is  a  multi-layered  task.  Variables  here  include, 
besides  the  different  tests,  sample  size  and  alternative  hypotheses.  That  is,  different  tests  will 
exhibit  different  powers  against  different  types  of  alternatives,  and  as  the  sample  size  varies.  It  is 
certainly  not  our  intent  here  to  engage  in  a  serious  study  and  ranking  of  the  many  possible  tests 
for  uniformity  of  spectrum.  Rather  we  shall  just  report,  as  a  rough  guide,  results  from  two  such 
tests  in  the  case  of  Methods  2  and  3. 
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We  first  consider  the  test  based  on  the  least  eigenvalue  which,  under  the  null  hypothesis  of 
uniform  spectrum,  has  the  distribution  given  by  Equation  (3.3).  From  each  of  the  two  methods 
for  generating  TRCMs,  we  obtain  a  sample  of  size  n  =  100,  for  TV  =  5,  10,  and  20,  compute 
the  spectrum,  and  then  the  empirical  distribution  function  of  the  least  eigenvalue.  This  is  then 
compared  with  the  theoretical  distribution  in  a  Kolmogorov  Smirnov  one-sample  test. 

In  the  case  of  TRCMs  generated  by  Method  2,  this  test  decisively  rejected  uniformity  of 
spectrum.  The  K-S  test  statistic  is  highly  significant,  as  shown  by  the  values  on  the  top  line 
of  Table  3-1.  The  problem  is  that  the  least  eigenvalue  is  too  large,  as  can  be  seen  in  the  case 
TV  =  10  from  Figure  3-1.  It  displays  the  empirical  and  theoretical  distributions,  and  also  indicates 
that  the  average  least  eigenvalue  =  0.2596,  with  a  sample  standard  deviation  of  0.1383.  Under 
the  null  hypothesis,  the  expected  value  of  the  least  eigenvalue  is  1/TV  =  0.1,  here,  with  vari¬ 
ance  =  (TV  —  1)/TV^(TV  -f  1)  =  0.00818,  here.  Hence  the  sample  average  should  be  approximately 
TV(0.1, 0.0000818)  distributed,  which  is  clearly  incompatible  with  the  data. 


TABLE  3-1 

K-S  Test  Statistics  for  Minimum 

Eigenvalue  of  TRCM 

Method  2 

N  =  5 

N  =  10 

N  =  20 

5.42 

5.54 

7.75 

Method  3 

1.91 

1.96 

2.05 

The  case  of  TRCMs  generated  by  Method  3  is  more  ambiguous,  partly  because  of  the  choice 
of  the  parameters  L  and  cr.  Certainly  there  are  ‘bad’  choices  of  these  parameters.  We  have  only 
mildly  tried  to  optimize  this  choice:  the  data  reported  in  the  lower  half  of  Table  3-1  is  based  on 
the  case  L  =  3TV  and  a  =  7TV.  We  see  that  the  K-S  statistic  is  still  significant  (values  >  1.63  are 
significant  at  the  1%  level).  But  the  sample  mean  is  only  about  2(7  away  from  its  mean  under 
the  null  hypothesis.  Overall,  this  case  presents  a  much  less  emphatic  rejection  of  uniformity  that 
associated  with  Method  2.  See  Figure  3-2  for  comparison. 

Our  second  test  involves  the  use  of  Neyman’s  statistic  TV2.  The  theory  of  this  test  is  available 
in  several  sources,  for  example  [19].  In  general,  there  is  a  sequence  {Nk  :  k  =  1,2,..}  of  statistics 
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Figure  3-L  Frequency  distribution  of  minimum  eigen vaiues:  Random  Toeplitz 
Method  2. 


that  can  be  used  to  test  the  hypothesis  of  uniformity  on  [0, 1]  against  an  alternative  density  of  the 
form 


f{x)  =  cexp 


1  + 

j=i 


where  the  Ij  are  Legendre  polynomials,  and  c  is  a  normalizing  constant  whose  precise  value  depends 
on  {bj}.  The  null  hypothesis  is  Hq  :  bj  =  0,  for  all  j,  or,  equivalently, 

i=i 

The  Neyman  statistic  Nk  =  uj  +  . . ,  +  where 

Vj  = 
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and  (x,  :  ^  =  1, . . . ,  n)  is  the  sample.  It  follows  that  when  /c  =  2,  A^2  is  a  combination  of  the  sample 
mean  and  variance  and  thus  has  an  intuitive  appeal.  It  is  known  that  each  is  asymptotically 
distributed  as  that  tests  b2ised  on  which  result  in  rejection  of  the  null  hypothesis 

for  large  values  of  Nk<,  are  consistent  and  asymptotically  unbiased.  In  practice,  the  x^(s)  approx¬ 
imation  to  N2  is  quite  accurate  for  n  >  20;  for  small  n,  upper  tail  percentage  point  tables  are 
available. 


NUMBER  OF  DATA 
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Figure  3-2.  Frequency  distribution  of  minimum  eigenvalues:  ilandom  Toeplitz 
Method  3. 


We  computed  the  Neyman  statistics  N2  for  the  same  batches  of  TRCMs  as  were  tested  using 
the  least  eigenvalue  statistic.  The  results  are  summarized  in  Table  3-2  below.  For  each  of  Methods  2 
and  3,  and  each  dimension  =  5,  10,  and  20,  we  give  the  average  and  standard  deviation  of  the 
100  N2's  computed  from  each  TRCM,  and  also  the  empirical  percentage  of  the  sample  whose  N2  - 
value  exceeded  the  .01  significance  level.  This  level,  for  the  cases  =  5,  10,  and  20  is,  respectively, 
9.64,  9.26,  and  9.23. 

The  main  inference  to  be  drawn  from  this  data  is  that  as  the  matrix  size  N  increases  beyond 
10,  we  observe  ever  more  pronounced  deviations  from  uniformity,  so  we  would  have  to  reject  the 
hypothesis  of  uniformity  for  TRCMs  generated  by  either  of  these  methods  for  N  >  10. 
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TABLE  3-2 

Neyman  N2  Statistic  Summary  for  TRCM 

N  =  5 

N  =  10 

N  =  20 

Ave. 

1.62 

4.16 

9.10 

S.D. 

0.86 

1.94 

3.30 

Method  2 

Sig. 

0% 

3% 

39% 

0.01 

Ave. 

2.62 

5.53 

11.50 

S.D. 

0.99 

1.69 

2.52 

Method  3 

Sig. 

0% 

3% 

83% 

0.01 

Our  conclusion,  based  on  both  tests,  is  that  for  N  <  10,  Method  3,  with  proper  choice  of  its 
parameters  L  and  cr,  can  be  said  to  produce  TRCMs  with  uniform  spectrum.  It  may  be  possible,  for 
larger  N,  to  adjust  these  parameters  in  Method  3,  so  as  to  pass  these  tests  for  uniform  spectrum. 
Failing  that,  and  insofar  as  a  random  spectrum  is  deemed  desirable,  we  must  fall  back  on  the  as 
yet  untried  Method  5. 
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4.  LOOSE  ENDS  AND  OPEN  QUESTIONS 


In  this  final  section  we  collect  together  a  number  of  remarks  pertinent  to  the  foregoing  mate¬ 
rial,  and  unresolved  issues  that  appear  worthy  of  further  study.  These  latter  issues  are  certainly 
of  varying  degrees  of  difficulty,  running  as  they  do  from  further  simulation  work  to  rather  deep 
theoretical  issues. 


4.1  Approximation  of  Toeplitz  Correlation  Matrices 


Given  that  our  primary  interest  in  TRCMs  is  that  they  serve  as  the  covariance  matrices  of 
sections  of  stationary  stochastic  processes,  it  is  conceivable  that  simple  models  for  such  processes 
might  yield  an  adequate  supply  of  TRCMs.  Of  course,  the  terms  ‘simple’  and  ‘adequate’  must  be 
carefully  specified. 


Let  us  say  that  ‘simple’  means  an  autoregressive  process  of  some  order  m,  which  we  denote 
by  AR{m),  The  cases  m  <  2  are  of  particular  interest.  When  m  =  1,  the  corresponding  N  x  N 
correlation  matrices  Cp  are  indexed  by  the  real  parameter  p,  |p|  <  1,  and  the  eigenvalues  of  Cp  are 
known  [17]  to  be  of  the  form 


Afc  = 


1 


I  <k<N, 


1  —  2p(:osu)k  H-  p^  ’ 

where  the  uJk  must  be  obtained  numerically  from  a  nonlinear  equation.  When  m  =  2,  the  eigenvalues 
are  not  similarly  available  in  closed  form,  but  can  still  be  indexed  by  a  pair  of  real  parameters.  For 
instance,  correlation  matrices  of  the  form  given  by  Equation  (2.4),  with 


4>{t)  =  exp(—  oc  1^1)  cos27rj0^, 

can  occur.  (In  m  =  1  case,  the  cosine  term  would  be  absent.) 

By  ‘adequate’  we  might  mean  either  the  degree  of  density  or  coverage  of  the  class  of  all  A  x  A 
Toeplitz  correlation  matrices  by  those  corresponding  to  the  AR[m)  processes,  or  the  degree  of 
coverage  of  the  eigenvalue  simplex  Sat  by  the  spectra  of  these  matrices.  Can  either  of  these  degrees 
of  coverage  by  correlation  matrices  of  AR{7n)  be  quantified  as  a  function  of  m?  And  if  so,  is  it 
small  enough  to  permit  usage  of  these  special  matrices  for  the  intended  applications?  But  the  real 
problem  may  be  that,  even  if  the  previous  questions  can  be  answered  affirmatively  for  some  m  <  A, 
it  will  prove  to  be  overly  difficult  to  actually  generate  the  correlation  matrices  corresponding  to  the 
appropriate  AR{m)  processes,  relative  to  the  methods  of  Section  3, 

We  remark  that  rational  approximation  in  the  frequency  domain  may  be  a  plausible  approach 
to  the  first  question  above. 


4.2  Length  of  Random  Vector  in  Method  2 

In  this  method  the  key  step  was  a  random  draw  from  the  unit  sphere  in  R^ .  It  is  possible  to 
draw  instead  from  the  unit  sphere  in  a  larger  space  L  >  N ^  and  still  define  the  function  (f)  as 
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before.  One  might  ask  how  the  spectral  properties,  in  particular  the  spectral  distribution  of  the 
resulting  TRCM,  depend  on  this  parameter  L.  A  brief  look  at  this  problem  has  suggested  that  the 
distribution  of  eigenvalues  does  not  become  more  uniform,  contrary  to  both  hope  and  expectation. 


4.3  Random  Probability  Distributions 

This  probably  is  the  most  interesting  theoretical  concept  to  arise  in  the  present  study.  In 
the  context  of  Method  3  we  have  discussed  one  possible  approach  to  this  problem,  namely,  the 
use  of  the  Dirichlet  process.  We  noted  that  conventional  methods  as  well  as  our  ad  hoc  approach 
yielded  discrete  probability  measures.  Here  we  want  to  mention  another,  more  recent,  procedure 
for  selecting  a  probability  measure  on  R  at  random  that,  with  probability  one,  yields  absolutely 
continuous  measures  with  an  analytic  density.  So,  in  a  sense,  we  are  at  the  opposite  extreme  from 
the  Dirichlet  distributions. 

This  procedure  is  due  to  Chen  and  Rubin  [4,  5],  and  is  based  on  the  use  of  an  orthonormal 
basis  for  L^{R).  (Other  measure  spaces  besides  (R,  Leb.)  could  be  used,  but  as  these  are  taken 
more  abstract,  some  of  the  special  properties  known  for  the  resulting  densities,  such  as  analyticity, 
existence  of  moments,  etc.  may  be  lost.)  If  {<^n}  is  such  a  basis,  and  b  =  {bn}  a  sequence  of  scalars 
such  that  Y^b^  =  1  (notation:  beS{P))^  then  the  function: 

fix)  =  (4.1) 

n 


is  easily  seen  to  be  a  probability  density  function  (pdf)  on  R.  It  follows  that  any  procedure  for 
selecting  sequences  b  randomly  from  the  sphere  S{P)  will  lead  to  a  random  pdf  via  Equation  (4.1). 

Of  course  there  is  a  difficulty  here  in  that  the  Hilbert  space  is  of  infinite  dimension  and 
hence  its  spheres  are  noncompact.  Further,  there  is  no  rotationally  invariant  probability  measure 
on  5(Z^),  although  there  are  finitely  additive  measures  on  the  algebra  of  cylinder  sets  of  which 
are  rotationally  invariant,  for  instance,  the  canonical  Gauss  measure.  But  in  order  that  a  cylinder 
set  measure  can  be  extended  to  be  a  countable  additive  measure  on  a  Borel  set  of  its  covariance 
operator  must  be  nuclear  [1].  This  means  that  its  spectrum  is  composed  of  eigenvalues  converging 
rapidly  enough  to  zero  so  as  to  form  a  summable  sequence.  Hence  the  associated  measure  is  highly 
anisotropic. 

Actually,  there  is  a  second  difficulty  here  besides  the  conceptual  one  associated  with  the  idea  of 
a  uniform  distribution  on  the  sphere  5(Z^),  and  that  concerns  the  necessity  to  approximate  infinite 
dimensional  vectors  by  finite  dimensional  ones.  That  is,  in  practice  the  sum  in  Equation  (4.1) 
must  be  truncated.  Fortunately,  sums  of  generally  not  unreasonable  length  suffice  to  approximate 
‘most’  pdfs,  as  shown  in  [4].  There  remains  the  task  of  numerical  simulation  of  random  draws  from 
S{1^)  to  achieve  such  approximations.  One  scheme  suggested  in  [4]  yields  an  expected  dimension 
of  about  50.  One  can  imagine  also  other  schemes  for  doing  the  simulations. 

The  point  to  be  made  here  is  that  we  have  here  a  new  method  for  producing  random  pdfs 
which  in  turn  will  lead  via  Method  3  to  TRCMs,  and  we  do  not  have  any  information  at  this  time 
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as  to  how  the  spectral  properties  of  these  matrices  behave.  Note,  however,  a  practical  drawback 
to  the  use  of  these  densities  for  generating  TRCMs  in  that  much  more  computational  effort  must 
be  expended  in  order  to  compute  their  Fourier  transforms,  in  comparison  wdth  the  case  of  discrete 
probability  measures. 

A  final  comment  is  to  recall  that  orthonormal  series  have  been  long  used  as  one  of  several 
popular  approaches  to  nonparametric  probability  density  estimation.  That  is,  an  estimate  of  an 
unknown  pdf  /  is  constructed  by  using  a  given  sample  to  form  estimates  of  a  finite  number  of 
Fourier  coefficients  of  /  with  respect  to  the  given  basis.  This  method  dates  from  work  of  N. 
Chencov,  1962,  with  further  contributions  by  S.  Schw^artz,  G.  Watson,  R.  Kronmal  and  M.  Tarter, 
and  D.  Bosq,  inter  alia.  Possibly  familiarity  with  this  approach  to  density  estimation  served  as 
some  motivation  for  the  Chen-Rubin  approach  to  density  construction.  We  note  that  while  their 
densities  are  nonnegative  by  definition,  as  they  should  be,  the  conventional  pdf  estimators  can  well 
lack  this  property  (particularly  for  small  data  samples),  a  well-known  shortcoming  of  the  orthogonal 
series  method. 


4.4  Other  Issues  in  Method  3 

There  are  several  parameters  that  can  be  varied:  the  integer  L  (number  of  point  masses),  the 
variance  parameter  a  (spread  of  the  points  about  the  orgin),  and  the  points  {i*}  themselves.  We 
have  used  tj  =  j  in  Equation  (2.5),  but  other  choices  are  possible,  and  might  lead  to  improved 
(more  uniform)  spectral  distributions.  Also,  the  set  {tj}  might  be  randomized. 

There  seems  to  be  a  distinct  lack  of  theoretical  results  about  these  specific  questions,  although 
there  is  certainly  no  lack  of  theory  of  symmetric  Toeplitz  matrices.  We  have  in  mind  especially 
results  about  the  asymptotic  behavior  of  eigenvalues.  Such  theorems  might  help  shed  light  on  the 
spectral  behavior  of  large  TRCMs  generated  by  either  of  Methods  2  or  3. 

If  the  urn  scheme  of  Equation  (2.7)  is  used  to  build  up  discrete  Dirichlet-distributed  probability 
mecLSures,  how  do  the  resulting  TRCMs  behave? 


4.5  Effect  of  the  Toeplitz  Projection 

Recall  that  the  projection  operation  C  — ►  Ct  described  in  Method  4  could  not  be  guaranteed 
to  ‘work,’  in  the  sense  of  always  returning  a  (Toeplitz)  correlation  matrix  Ct  whenever  a  correlation 
matrix  C  was  input.  Nevertheless  it  ‘usually’  works.  This  phenomenon  seems  worthy  of  further 
understanding  and  quantification.  Is  it  possible  to  assign  a  probability  to  the  statement 

[{CeT{Nyi 

More  generally,  one  might  inquire  if  there  are  any  general  facts  about  changes  in  spectrum 
induced  by  projections  in  the  space  of  square  matrices  with  the  Hilbert-Schmidt  inner  product  of 
Equation  (2.8).  A  particular  question  is  what  the  Toeplitz  projection  does  to  a  class  of  ROMs  with 
uniform  spectrum,  at  least  in  those  cases  where  it  preserves  the  positivity. 


19 


Apropos  of  this  last  question,  there  arises  a  certain  practical  matter:  how  does  one  produce 
a  sequence  of  RCMs  with  a  uniform  spectrum?  Note  that  this  is  not  a  theoretical  difficulty  — 
one  begins  with  any  positive  definite  matrix  with  a  desired  spectrum  in  5^^  and  then  sequentially 
adjusts  the  diagonal  entries  via  a  sequence  of  Givens  rotations.  This  method  is  discussed  in  [11], 
with  alternate  references  provided,  and  it  is  noted  that  the  IMSL  subroutine  GGCOR  provides 
a  commercial  package  for  this  very  purpose  (this  subroutine  will  also  make  random  draws  from 
the  orthogonal  group  0{N)).  The  problem  is  that,  as  the  subroutine  documentation  warns,  the 
spectrum  of  the  final  matrix  may  differ  from  the  initial  input  spectrum.  In  fact,  this  problem  seems 
to  be  very  real.  In  100  trials  each  for  =  3,  5,  8,  with  an  initial  spectrum  drawn  randomly  from 
the  uniform  distribution  [7(5^'),  we  observed  an  empirical  failure  rate  of  0.3%  (resp.,  7%,  18%). 
That  is,  these  percentages  of  the  trials  resulted  in  nonpositive  definite  output  matrices. 


4.6  Toeplitz  RCMs  with  Random  Spectrum 

We  have  seen  that,  especially  for  matrices  of  order  N  >  10,  we  do  not  have  as  yet  a  method 
of  producing  Toeplitz  RCMs  with  a  random  (uniform)  spectrum  other  than  Method  5,  with  which 
we  have  no  direct  experience.  We  have  mentioned  the  baisic  questions  of  existence  and  convergence 
associated  with  this  method,  and  the  probable  difficulties  with  implementation  and  round-off. 

Another  basic  question  which  needs  to  be  recorded  in  this  context  is  whether  we  should  expect 
the  existence  of  TRCMs  with  random  spectrum.  Our  numerical  work  leads  us  to  accept  this 
hypothesis  for  small  orders  {N  <  10),  but  is  otherwise  inconclusive.  But,  in  any  event,  there  is  no 
existence  theorem,  and  the  connection  between  such  a  result,  if  true,  and  the  inverse  eigenvalue 
problem  for  Toeplitz  correlation  matrices  should  be  explored.  Certainly,  if  the  inverse  eigenvalue 
problem  is  solvable,  and  the  convergence  of  Method  5,  or  some  other,  is  validated,  then  we  will 
have  a  source  of  TRCMs  with  random  spectrum.  And,  in  the  converse  direction,  if  TRCMs  with 
random  spectrum,  can  be  constructed,  then,  with  arbitrarily  high  probability,  we  can  find  a  Toeplitz 
correlation  matrix  with  a  given  spectrum.  (Sketch  of  argument  ...  Let  0  <  Pn  T  fn  be  a 

decreasing  sequence  of  open  sets  in  the  simplex  Sat  with  nVn  =  {Aq}.  Then  with  probability  >  pn 
there  is  a  Toeplitz  correlation  matrix  Tn  with  a{Tn)sVn-  By  compactness  of  the  set  of  N  x  N 
correlation  matrices,  there  is  a  convergent  subsequence  of  {Tn}  with  limit  T,  and  by  standard 
results  on  spectral  perturbations,  (t{T)  =  {Aq}.  The  probability  that  the  first  m  of  the  Tn  exists  is 
Pi  •  •  -  Pm  and,  if  the  increase  of  pn  to  1  is  sufficiently  rapid,  this  partial  product  will  coverage  to  a 
number  arbitrarily  close  to  1.) 


4.7  Further  Simulation  and  Power  Studies 

Failing  theoretical  advances,  there  is  always  the  possibility  of  further  numerical  work  along 
the  lines  reported  in  Section  3.  This  might  take  the  form  of  either  additional  tests  for  uniformity 
or,  more  interestingly,  studies  comparing  the  power  of  various  tests  for  uniformity  relative  to  a 
particular  alternative.  In  particulair,  is  there  a  uniformly  most  powerful  test  for  uniformity? 
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As  an  example  we  looked  at  a  statistic  of  Lohrding  [14]  which  was  shown  to  have  favorable 
power  against  various  alternatives  with  respect  to  other  popular  test  statistics  such  as  Kolmogorov- 
Smirnov  and  Cramer-von  Mises.  In  the  case  where  we  are  testing  an  ordered  sample  {xi, . . . , } 
for  uniformity  on  [0, 1],  the  test  statistic  in  question  is 


T  —  {d\  +  •  ’  •  +  —  1), 


where 


dk  —  \^k  ^k\/^ki 

here  and  are  the  mean  and  variance  of  the  kth  order  statistic  from  a  sample  of  size  TV  —  1 
from  U[0,  Ij.  Use  of  this  statistic  for  TRCMs  generated  from  Methods  2  and  3  bcisically  reconfirmed 
results  reported  in  Section  3  for  the  other  two  tests,  by  leading  to  rejection  of  the  null  hypothesis 
for  TV  >  10. 
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5.  SUMMARY 


The  theory  and  practice  of  random  Toeplitz  correlation  matrices  is  more  subtle  than  that  of 
random  correlation  matrices.  This  is  to  be  expected  given  the  corresponding  greater  mathemat¬ 
ical  difficulties  with  Toeplitz  matrices  and  operators  vis-a-vis  Hermitian  matrices  and  operators. 
Starting  from  this  background,  we  prepared  this  report  with  two  goals  in  mind:  to  raise  several  the¬ 
oretical  issues  concerning  the  spectral  behavior  of  Toeplitz  RCMs,  and  to  present  some  numerical 
experience  with  certain  special  cases. 

In  fact,  a  variety  of  unresolved  theoretical  questions  arose  in  the  course  of  this  work,  and  these 
are  collected  together  in  Section  4.  The  most  interesting  of  these  appear  to  be  (1)  the  spectral 
behavior  of  Toeplitz  RCMs  generated  from  the  characteristic  functions  of  various  random  proba¬ 
bility  measures;  (2)  the  connection  between  the  inverse  eigenvalue  problem  for  Toeplitz  correlation 
matrices  and  the  concept  of  Toeplitz  RCMs  with  random  (uniform)  spectrum;  and  (3)  the  behavior 
of  the  matrix  spectrum  under  the  Toeplitz  projection  operator. 

The  numerical  results  shed  varying  degrees  of  light  on  some  of  the  theoretical  questions,  as  we 
have  tried  to  indicate.  Insofar  as  the  matter  of  computer  generation  of  Toeplitz  RCMs  is  concerned, 
we  are  left  with  Method  3  (based  on  the  characteristic  functions  of  random  discrete  probability 
measures)  as  the  method  of  choice.  This  method  is  not  quite  as  efficient  as  Method  2  (based  on  the 
autocorrelation  function  of  a  random  discrete  moving  average  process),  but  seems  to  yield  more 
uniformly  distributed  spectra.  But  the  tests  of  both  methods  for  uniformity  fail  for  moderately 
large  matrix  sizes  {N  >  10),  and  so  that  area  remains  open  too. 
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